Step selection simulations
Generating simulations from a fitted (temporally dynamic) step selection function
Scott Forrest
2024-02-29
In this script we will simulate trajectories from a temporally dynamic step selection function (SSF) fitted to GPS data. The model fitting is done in a separate script, as well as the construction of the temporally varying coefficients from the harmonic terms.
Load packages
options(scipen=999)
library(tidyverse)
packages <- c("amt", "lubridate", "terra", "tictoc",
"beepr", "matrixStats", "Rfast")
walk(packages, require, character.only = T)Load a table of temporally varying coefficients (reconstructed from the harmonic terms) from the model fitting script. This has coefficients at 0.1 hour increments, but we are simulating steps at 1 hourly intervals, so we will subset to retain coefficients at 1 hourly intervals.
We are using the coefficients from the model with three pairs of harmonics. To use coefficients from a different model, the only change in the script should be the table of coefficients.
These coefficients should be on the natural scale, as they will be multiplied by the habitat covariates in their natural scales. For more info on this see the model fitting script.
We will plot the temporally varying coefficients to check whether everything looks as it should.
#Load coefficient estimates from the model fitting script
hourly_coefs <- read_csv("outputs/TwoStep_3pDaily_coefs_dry_2024-02-20.csv")## Rows: 240 Columns: 15
## ── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (15): hour, ndvi, ndvi_2, canopy, canopy_2, slope, herby, memory, memory_2, sl, log_sl, cos_ta, shape, scale, kappa
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(hourly_coefs)## # A tibble: 6 × 15
## hour ndvi ndvi_2 canopy canopy_2 slope herby memory memory_2 sl log_sl cos_ta shape scale kappa
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 8.42 -14.4 1.84 -3.52 -0.223 0.149 21297718. -3.12e13 -0.00169 0.0632 -0.173 0.501 281. 0.0121
## 2 0.1 8.53 -14.5 1.85 -3.53 -0.227 0.147 20905425. -3.00e13 -0.00178 0.0680 -0.176 0.506 274. 0.00881
## 3 0.2 8.62 -14.6 1.85 -3.54 -0.231 0.145 20514573. -2.88e13 -0.00188 0.0725 -0.179 0.511 266. 0.00536
## 4 0.3 8.71 -14.7 1.85 -3.55 -0.235 0.143 20125640. -2.76e13 -0.00198 0.0766 -0.183 0.515 260. 0.00170
## 5 0.4 8.79 -14.8 1.85 -3.56 -0.239 0.140 19739083. -2.65e13 -0.00207 0.0802 -0.187 0.518 254. -0.00222
## 6 0.5 8.86 -14.8 1.85 -3.57 -0.242 0.137 19355341. -2.53e13 -0.00216 0.0834 -0.191 0.522 248. -0.00644# lengthen data frame for ggplot
hourly_coefs_long <- hourly_coefs %>% pivot_longer(cols = !hour)
# uncomment variables to also plot them, but some of these (particularly memory)
# are on very different scales (as the previous space use density is tiny)
hourly_coefs_long %>%
filter(!name %in% c("shape", "scale", "kappa", "memory", "memory_2")) %>%
ggplot(aes(x = hour, y = value, colour = factor(name))) +
geom_line() +
scale_colour_viridis_d("Covariate") +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_x_continuous("Hour", breaks = seq(0,24,2)) +
theme_classic()
Subset to retain coefficients at 1 hourly intervals. We’ll also set hour 0 to be hour 24 (as they are the same thing) so we can index from 1.
Again we will plot the temporally varying coefficients to check whether the subsetting worked correctly.
# keep only the integer hours using the modulo operator, which returns the remainder of the division of the first argument by the second
hourly_coefs <- hourly_coefs %>% filter(hourly_coefs$hour %% 1 == 0) %>%
mutate(hour = ifelse(hour == 0, 24, hour))
head(hourly_coefs)## # A tibble: 6 × 15
## hour ndvi ndvi_2 canopy canopy_2 slope herby memory memory_2 sl log_sl cos_ta shape scale kappa
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 24 8.42 -14.4 1.84 -3.52 -0.223 0.149 21297718. -3.12e13 -0.00169 0.0632 -0.173 0.501 281. 0.0121
## 2 1 9.06 -15.1 1.86 -3.66 -0.260 0.121 17493049. -2.00e13 -0.00260 0.0910 -0.218 0.529 224. -0.0336
## 3 2 8.60 -14.8 1.85 -3.83 -0.281 0.0821 14195313. -1.16e13 -0.00321 0.0613 -0.309 0.499 197. -0.124
## 4 3 7.18 -13.6 1.77 -3.86 -0.275 0.0475 11769173. -7.39e12 -0.00334 -0.0115 -0.428 0.427 192. -0.243
## 5 4 5.43 -12.1 1.61 -3.61 -0.240 0.0252 10666304. -7.72e12 -0.00289 -0.0779 -0.504 0.360 210. -0.319
## 6 5 4.28 -11.2 1.45 -3.06 -0.188 0.0103 11431084. -1.28e13 -0.00193 -0.0826 -0.453 0.356 263. -0.268# lengthen data frame for ggplot
hourly_coefs_long <- hourly_coefs %>% pivot_longer(cols = !hour)
# uncomment variables to also plot them, but some of these (particularly memory)
# are on very different scales (as the previous space use density is tiny)
hourly_coefs_long %>%
filter(!name %in% c("shape", "scale", "kappa", "memory", "memory_2")) %>%
ggplot(aes(x = hour, y = value, colour = factor(name))) +
geom_line() +
scale_colour_viridis_d("Covariate") +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_x_continuous("Hour", breaks = seq(0,24,2)) +
theme_classic()
Load the memory parameters (to determine previous space use density)
from the Memory_parameter_estimation_ scripts. For the
simulations in the paper we used the
Memory_parameter_estimation_loc-subset.Rmd script. The file
that we are reading in contains the mean of the KDE
bandwidth, \(\bar{\sigma}\) (Gaussian
kernel standard deviation), and the temporal exponential decay parameter
for the optimisation over all individuals, \(\gamma_{pop}\).
The Memory_parameter_estimation_ scripts also outputs a
csv that has individual-level KDE bandwidths, \(\sigma_i\), and temporal exponential decay
parameters, \(\gamma_i\) where \(i\) is an individual buffalo.
The parameters that relate to the model that was fitted should be used - if the individual-level memory parameters were used to estimate the previous space use density and then individual SSFs were fitted, then the individual-level memory parameters should be used in the simulations. If the population-level memory parameters were used to estimate the previous space use density and then a population-level SSF was fitted, then the population-level memory parameters should be used in the simulations.
memory_params_ALLoptim <- read_csv("outputs/memory_params_ALLoptim_2024-02-05.csv")## New names:
## Rows: 1 Columns: 3
## ── Column specification
## ──────────────────────────────────────────────────────────────────────────────────────────────── Delimiter: "," dbl
## (3): ...1, exp_gamma_all_optim, mean_kde_sd
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ Specify the column types or set
## `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`memory_params_ALLoptim## # A tibble: 1 × 3
## ...1 exp_gamma_all_optim mean_kde_sd
## <dbl> <dbl> <dbl>
## 1 1 0.00737 572.Import the buffalo starting locations (to start the simulations from). This is simple enough to generate by importing data and selecting the first location, although we import the file here so we don’t have to import a large dataset.
start_buffalo_df <- read_csv("outputs/buffalo_starting_locations_13inds_2024-02-06.csv")## Rows: 13 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): id, x, y
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(start_buffalo_df)## # A tibble: 6 × 3
## id x y
## <dbl> <dbl> <dbl>
## 1 2005 41968. -1435673.
## 2 2014 33615. -1440916.
## 3 2018 30093. -1423073.
## 4 2021 25964. -1437033.
## 5 2022 27928. -1421152.
## 6 2024 8227. -1416224.Import environmental raster layers that the SSF model was fitted with.
# unscaled rasters
ndvi_stack <- rast("mapping/cropped rasters/ndvi_GEE_projected_watermask20230207.tif")
canopy <- rast("mapping/cropped rasters/canopy_cover.tif")
herby <- rast("mapping/cropped rasters/veg_herby.tif")
slope <- rast("mapping/cropped rasters/slope_raster.tif")Predictions
The external selection of covariates is denoted by
\[ \omega(\mathbf{X}(s_t); \boldsymbol{\beta}(\tau; \boldsymbol{\alpha})) = \exp(\beta_{1}(\tau; \boldsymbol{\alpha}_1) X_1(s_t) + \cdots + \beta_{n}(\tau; \boldsymbol{\alpha}_n) X_n(s_t)), \]
where \(\mathbf{X}(s_t)\) is a vector of covariate values at location \(s_t\), \(\boldsymbol{\beta}(\tau; \boldsymbol{\alpha})\) is the vector of temporally varying coefficients, which is a function of time \(\tau\), and \(\boldsymbol{\alpha}\), which is the vector of parameters defined by the harmonic terms.
To aid the computation of hte simulations, we can precompute \(\omega(\mathbf{X}(s_t); \boldsymbol{\beta}(\tau; \boldsymbol{\alpha}))\) prior to running the simulations.
In the dataframe of temporally varying coefficients, for each covariate we have reconstructed \(\beta_{i}(\tau; \boldsymbol{\alpha}_i)\) and discretised for each hour of the day, resulting in \(\beta_{i,\tau}\) for \(i = 1, \ldots, n\) where \(n\) is the number of covariates and \(\tau = 0, \ldots, 23\).
Given these, we can solve \(\omega(\mathbf{X}(s_t); \boldsymbol{\beta}(\tau; \boldsymbol{\alpha}))\) for every hour of the day. This will result in a spatial map of the selection of covariates for each hour of the day.
Then, when we do our step selection simulations, we can just subset these maps by the current hour of the day, and extract the values of \(\omega(\mathbf{X}(s_t); \boldsymbol{\beta}(\tau; \boldsymbol{\alpha}))\) for each proposed step location, rather than solving \(\omega(\mathbf{X}(s_t); \boldsymbol{\beta}(\tau; \boldsymbol{\alpha}))\) for every step location.
Calculating \(\omega(\mathbf{X}(s_t); \boldsymbol{\beta}(\tau; \boldsymbol{\alpha}))\) has been called the `naive’ prediction approach by Signer, Fieberg, and Avgar (2017) as it doesn’t account for the movement dynamics, which we will incorporate through the simulations.
Naive prediction approach
As we fitted either wet or dry season models to data that covered 2018 and 2019, we will use the mean of the NDVI for the parts of the 2018 and 2019 dry season that the data covered. We will also square NDVI and canopy cover now, and then it’s just a multiplication with the temporally varying coefficient within the function.
ndvi_2018_dry <- ndvi_stack[[8:10]]
ndvi_2019_dry <- ndvi_stack[[19:21]]
ndvi <- terra::mean(c(ndvi_2018_dry, ndvi_2019_dry))
names(ndvi) <- "NDVI"
plot(ndvi)
# for plotting with ggplot after generating simulations
ndvi_df <- as.data.frame(ndvi, xy = TRUE)
ndvi_df$NDVI_discrete <- cut(ndvi_df$NDVI, breaks=9, dig.lab = 2)
# NDVI squared
ndvi_sq <- ndvi ^ 2
# rescale canopy cover to 0 - 1 (as this was what the model was fitted to)
canopy01 <- canopy/100
# canopy cover squared
canopy01_sq <- canopy01 ^ 2Generating the ‘naive’ habitat selection layers for stochastic simulations
# creaty empty objects to store the layers
naive_pred_stack <- c()
tic("Naive predictions")
# for each hour of the day
for(hour_no in 1:24) {
# create a raster stack of the covariates (including the squared terms)
resources <- c(ndvi,
ndvi_sq,
canopy01,
canopy01_sq,
slope,
herby
)
# ndvi
# using the linear term
ndvi_lin <- resources[[1]] *
hourly_coefs$ndvi[[which(hourly_coefs$hour == hour_no)]]
# using the quadratic term
ndvi_quad <- resources[[2]] *
hourly_coefs$ndvi_2[[which(hourly_coefs$hour == hour_no)]]
# combining
ndvi_pred <- ndvi_lin + ndvi_quad
# canopy cover
# using the linear term
canopy_lin <- resources[[3]] *
hourly_coefs$canopy[[which(hourly_coefs$hour == hour_no)]]
# using the quadratic term
canopy_quad <- resources[[4]] *
hourly_coefs$canopy_2[[which(hourly_coefs$hour == hour_no)]]
# combining
canopy_pred <- canopy_lin + canopy_quad
# herby
slope_lin <- resources[[5]]
slope_pred <- slope_lin *
hourly_coefs$slope[[which(hourly_coefs$hour == hour_no)]]
# herby
herby_lin <- resources[[6]]
herby_pred <- herby_lin *
hourly_coefs$herby[[which(hourly_coefs$hour == hour_no)]]
# combining all covariates (but not exponentiating yet)
naive_pred <- ndvi_pred + canopy_pred + slope_pred + herby_pred
# adding to the list of rasters for each hour
naive_pred_stack <- c(naive_pred_stack, naive_pred)
}
# plot some example naive predictions
plot(naive_pred_stack[[3]], main = "Naive predictions - hour 3")
plot(naive_pred_stack[[12]], main = "Naive predictions - hour 12")
# if we want to save the stack to reimport later so we don't have to recompute
filename <- paste0("mapping/CLR_TS_2pDailyMemALL_naive_unnorm_hourly_",
Sys.Date(), ".tif")
# filename
terra::writeRaster(rast(naive_pred_stack), filename, overwrite = TRUE)##
|---------|---------|---------|---------|
=========================================
toc()## Naive predictions: 91.49 sec elapsedbeep(sound = 2)Cropping the predictions and setting the origin
If using a ‘wrapped’ or ‘torus’ boundary, then the origin should be
set to (0,0) as we’ll use the modulo operator (which is %%
in R) to do the wrapping. We will set it back to it’s original
coordinates after running the simulations. To be consistent we will just
always set the origin to (0,0).
We can also crop layers at this point if we want a smaller extent.
Just change the xmin, xmax, ymin,
ymax values. Here we won’t crop, so we just pull the extent
values from the layer itself, but we do use these values (actually just
xmin and ymin) to set the origin.
# convert the list of naive prediction rasters to a raster stack
naive_predictions <- rast(naive_pred_stack)
# set the time of the rasters - here we are only interested in the hour,
# so we just set an arbitrary date. Although if predicting over seasons for
# instance, then set the date with the correct increment
terra::time(naive_predictions) <- ymd_hm("2018-08-01 00:00",
tz = "Australia/Queensland") + hours(1:24)
terra::time(naive_predictions)## [1] "2018-08-01 01:00:00 AEST" "2018-08-01 02:00:00 AEST" "2018-08-01 03:00:00 AEST" "2018-08-01 04:00:00 AEST"
## [5] "2018-08-01 05:00:00 AEST" "2018-08-01 06:00:00 AEST" "2018-08-01 07:00:00 AEST" "2018-08-01 08:00:00 AEST"
## [9] "2018-08-01 09:00:00 AEST" "2018-08-01 10:00:00 AEST" "2018-08-01 11:00:00 AEST" "2018-08-01 12:00:00 AEST"
## [13] "2018-08-01 13:00:00 AEST" "2018-08-01 14:00:00 AEST" "2018-08-01 15:00:00 AEST" "2018-08-01 16:00:00 AEST"
## [17] "2018-08-01 17:00:00 AEST" "2018-08-01 18:00:00 AEST" "2018-08-01 19:00:00 AEST" "2018-08-01 20:00:00 AEST"
## [21] "2018-08-01 21:00:00 AEST" "2018-08-01 22:00:00 AEST" "2018-08-01 23:00:00 AEST" "2018-08-02 00:00:00 AEST"hour(time(naive_predictions))## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 0# full extent
xmin <- ext(naive_predictions)[1]
xmax <- ext(naive_predictions)[2]
ymin <- ext(naive_predictions)[3]
ymax <- ext(naive_predictions)[4]
crop_extent <- ext(xmin, xmax, ymin, ymax)
# if we wanted to set the extent just around the buffalo data (with a buffer)
# buffer = 5000 # in metres
# xmin <- round(min(buffalo_data$x2), digits = -2) - buffer
# xmax <- round(max(buffalo_data$x2), digits = -2) + buffer
# ymin <- round(min(buffalo_data$y2), digits = -2) - buffer
# ymax <- round(max(buffalo_data$y2), digits = -2) + buffer
# crop_extent <- ext(xmin, xmax, ymin, ymax)
# to create a smaller extent centred around the centre point
# x_centrepoint <- (xmin+xmax)/2
# y_centrepoint <- (ymin+ymax)/2
# cell_size <- 25
# width <- 500 # number of cells
# height <- 500 # number of cells
# xmin <- x_centrepoint - (cell_size * width/2)
# xmax <- x_centrepoint + (cell_size * width/2)
# ymin <- y_centrepoint - (cell_size * height/2)
# ymax <- y_centrepoint + (cell_size * height/2)
# crop_extent <- ext(xmin, xmax, ymin, ymax)
# crop the naive predictions
naive_predictions_cropped <- terra::crop(naive_predictions, crop_extent)
# set the origin at (0,0)
ext(naive_predictions_cropped) <- c(xmin - xmin, xmax - xmin,
ymin - ymin, ymax - ymin)Setting up the simulation function
There are several parts to the memory process as it needs to ‘warm-up’ so there are enough locations to start using as the memory. This part of the trajectory is then discarded.
simulate_ssf_memory <- function(
n_steps, # final number of steps in the trajectory
n_ch, # number of proposed steps
coef, # table of temporally varying coefficient values
xy0, # starting location
starting_time, # time of the first step in POSIXct format
step_time_increment, # time increment between steps in minutes
resc, # resources - naive predictions
# boundary = "wrapped", # boundary type (either "wrapped" or "reflective")
memory_period, # number of locations to be included in the temporal decay
memory_delay, # number of locations to exclude prior to the current step
spatial_sd, # KDE bandwidth parameter for estimating previous space use
exp_gamma # temporal decay parameter for estimating previous space use
) {
tic()
# vectorisation of the random sampling leads to iterating over hourly coef values
sl <- rgamma(n_steps * n_ch,
shape = rep(coef$shape, each = n_ch),
scale = rep(coef$scale, each = n_ch))
# create a vector of the turning angle coefficients (over the full trajectory)
ta_coefs <- rep(coef$kappa, each = n_ch, length.out = n_steps)
# if positive, mean = 0 (which is rvonmises(mu = pi) - pi),
# indicating persistent steps in a similar direction,
# otherwise mean = pi (which is rvonmises(mu = 0) - pi),
# indicating reversal steps
ta_coefs_mu <- ifelse(ta_coefs > 0, pi, 0)
# sample turning angles with the different kappa and mean values
ta <- as.vector(mapply(Rfast::rvonmises,
n = n_ch,
m = ta_coefs_mu,
k = abs(ta_coefs))) - pi
# setup the simulation
steps <- rep(1:n_steps, each = n_ch)
step_time <- starting_time + minutes((0:(n_steps-1))*step_time_increment)
# set the starting location and (random) angle
x_0 <- xy0[[1]]
y_0 <- xy0[[2]]
angle_0 <- runif(1, min = -pi, max = pi)
# create vectors to store the results
x <- rep(NA, n_steps)
y <- rep(NA, n_steps)
hour <- rep(NA, n_steps)
step_length <- rep(NA, n_steps)
angle <- rep(NA, n_steps)
bearing <- rep(NA, n_steps)
hab_p <- rep(NA, n_steps)
mem_p <- rep(NA, n_steps)
# set initial values
x[1] <- x_0
y[1] <- y_0
hour[1] <- 1
step_length[1] <- 0
angle[1] <- angle_0
bearing[1] <- 0
hab_p[1] <- 0
mem_p[1] <- 0
for (i in 2:n_steps) {
# if you want the option for the boundary to be wrapped or reflective,
# uncomment the following lines, and then specify in the function arguments
# if(boundary == "wrapped") {
# wrapped boundary
# adding the modulo operator %% ensures a wrapped/toroid landscape,
# but the origin must be at (0,0)
# bearing <- atan2(y[i - 2] - y[i - 1], x[i - 2] - x[i - 1])
#
# x_prop <- (x[i - 1] + sl[steps == i] * cos(bearing[i - 1] +
# ta[steps == i])) %% ext(resc)[2]
#
# y_prop <- (y[i - 1] + sl[steps == i] * sin(bearing[i - 1] +
# ta[steps == i])) %% ext(resc)[4]
#
# angle_prop <- ta[steps == i]
# bearing_prop <- atan2(y_prop - y[i - 1], x_prop - x[i - 1])
# } else if(boundary == "reflective") {
#
# # reflective boundary
# proposed step in the x-direction
x_prop <- (x[i - 1] + sl[steps == i] * cos(bearing[i - 1] + ta[steps == i]))
# proposed step in the y-direction
y_prop <- (y[i - 1] + sl[steps == i] * sin(bearing[i - 1] + ta[steps == i]))
# proposed step length
sl_prop <- sl[steps == i]
# proposed turning angle
angle_prop <- ta[steps == i]
# proposed bearing (required for imposing correlated movement,
# as the turning angle of the following step is relative to the bearing)
bearing_prop <- atan2(y_prop - y[i - 1], x_prop - x[i - 1])
# hour of the proposed step, which could be day, month, etc for other time scales
hour_prop <- lubridate::hour(step_time[i])
#
# } else {
# print("Specify either 'reflective' (hard boundary) or
# 'wrapped' (toroidal) for the boundary")
# }
# sampling environmental suitability from the naive prediction rasters
# pull out the relevant raster based on the time and extract values from it
# these values will be exponentiated after added to the log-memory probability
p <- terra::extract(
resc[[which(lubridate::hour(terra::time(naive_predictions)) == hour_prop)]],
cbind(x_prop, y_prop))[,1]
### memory process ###
# the coefficients (including memory process) for hour span from 1 to 24
# rather than 0 to 23, so we'll update the hour_prop index to 24 if it's 0,
# as hour 0 and hour 24 are the same
hour_prop_mem <- ifelse(hour_prop == 0, 24, hour_prop)
# if the memory period is less than 2x the memory delay,
# then all steps are included and there is no memory delay
# (e.g. excluding the past 24 hours)
if(sum(!is.na(x)) < 2*memory_delay) {
memory_density <- mapply(function(x_prop, y_prop) {
# subset all previous steps
recent_x <- x[1:(i-1)]
recent_y <- y[1:(i-1)]
# calculate the distance and time between the current and previous steps
diff_x <- x_prop - recent_x
diff_y <- y_prop - recent_y
diff_time <- length(recent_x):1
# calculate the memory density
# the logSumExp function is used to avoid numerical underflow and overflow
# we calculate all components (i.e. the x and y spatial and temporal components)
# in the log-scale
# we normalise by the temporal decay component
return(logSumExp(dnorm(diff_x, mean = 0, sd = spatial_sd, log = TRUE) +
dnorm(diff_y, mean = 0, sd = spatial_sd, log = TRUE) +
(-exp_gamma*diff_time), na.rm = TRUE) -
logSumExp(-exp_gamma*diff_time, na.rm = TRUE))
}, x_prop, y_prop) # pass current locations in here down
# calculate the memory probability by multiplying the density
# (and it's square) by the coefficients
# for calculations it is
# better to use the log-scale above,
# so we exponentiate here to get it back to the natural scale (which was the
# scale of our covariate when fitting the model)
# this should not be confused with exponentiating the get the external
# selection probability, which we do below after adding the habitat suitability
# from the naive predictions
memory_p <- coef$memory[hour_prop_mem] * exp(memory_density) +
coef$memory_2[hour_prop_mem] * exp(memory_density)^2
# if the memory period is greater than 2x the memory delay,
# but less than the memory period, then all steps are considered,
# but there is a memory delay (e.g. excluding the past 24 hours)
} else if(sum(!is.na(x)) <= memory_period) {
memory_density <- mapply(function(x_prop, y_prop) {
# subset all previous steps, but exclude the memory delay
recent_x <- x[1:(i-memory_delay)]
recent_y <- y[1:(i-memory_delay)]
# calculate the distance and time between the current and previous steps
diff_x <- x_prop - recent_x
diff_y <- y_prop - recent_y
diff_time <- length(recent_x):1
# calculate the memory density
# the logSumExp function is used to avoid numerical underflow and overflow
# we calculate all components (i.e. the x and y spatial and temporal components)
# in the log-scale
# we normalise by the temporal decay component
return(logSumExp(dnorm(diff_x, mean = 0, sd = spatial_sd, log = TRUE) +
dnorm(diff_y, mean = 0, sd = spatial_sd, log = TRUE) +
(-exp_gamma*diff_time)
, na.rm = TRUE) -
logSumExp(-exp_gamma*diff_time, na.rm = TRUE))
}, x_prop, y_prop) # pass current locations in here down
# calculate the memory probability by multiplying the density
# (and it's square) by the coefficients
# for calculations it is
# better to use the log-scale above,
# so we exponentiate here to get it back to the natural scale (which was the
# scale of our covariate when fitting the model)
# this should not be confused with exponentiating the get the external
# selection probability, which we do below after adding the habitat suitability
# from the naive predictions
memory_p <- coef$memory[hour_prop_mem] * exp(memory_density) +
coef$memory_2[hour_prop_mem] * exp(memory_density)^2
# if the memory period is greater than the memory period,
#then only the memory period steps are considered (e.g. excluding the past 24 hours)
} else if(sum(!is.na(x)) > memory_period) {
memory_density <- mapply(function(x_prop, y_prop) {
# now we subset from the start of the memory period to the memory delay,
# i.e. a maximum of memory_period - memory_delay steps
recent_x <- x[(i-memory_period):(i-memory_delay)]
recent_y <- y[(i-memory_period):(i-memory_delay)]
# calculate the distance and time between the current and previous steps
diff_x <- x_prop - recent_x
diff_y <- y_prop - recent_y
diff_time <- length(recent_x):1
# calculate the memory density
# the logSumExp function is used to avoid numerical underflow and overflow
# we calculate all components (i.e. the x and y spatial and temporal components)
# in the log-scale
# we normalise by the temporal decay component
return(logSumExp(dnorm(diff_x, mean = 0, sd = spatial_sd, log = TRUE) +
dnorm(diff_y, mean = 0, sd = spatial_sd, log = TRUE) +
(-exp_gamma*diff_time)
, na.rm = TRUE) -
logSumExp(-exp_gamma*diff_time, na.rm = TRUE))
}, x_prop, y_prop) # pass current locations in here down
# calculate the memory probability by multiplying the density
# (and it's square) by the coefficients
# for calculations it is
# better to use the log-scale above,
# so we exponentiate here to get it back to the natural scale (which was the
# scale of our covariate when fitting the model)
# this should not be confused with exponentiating the get the external
# selection probability, which we do below after adding the habitat suitability
# from the naive predictions
memory_p <- coef$memory[hour_prop_mem] * exp(memory_density) +
coef$memory_2[hour_prop_mem] * exp(memory_density)^2
} else {
print("Out of specified memory conditions - check memory arguments")
}
# if the probability of either the habitat suitability probability or
# memory probability is NA (i.e. if falls outside the extent when using
# the reflective boundary, or in non-habitat such as water),
# then set it to a very low value (prior to exponenitating, e.g. -50)
p[is.na(p)] <- -50
memory_p[is.na(memory_p)] <- -50
# calculate the total probability of the proposed step and
# sample one of the locations based on this probability
w <- sample(n_ch, 1, prob = exp(memory_p + p))
# store all the information for the proposed step
x[i] <- x_prop[w]
y[i] <- y_prop[w]
hour[i] <- hour_prop
step_length[i] <- sl_prop[w]
angle[i] <- angle_prop[w]
bearing[i] <- bearing_prop[w]
hab_p[i] <- p[w]
mem_p[i] <- memory_p[w]
}
(toc())
return(data_frame(x = x,
y = y,
t = step_time,
hour = hour,
sl = step_length,
angle = angle,
bearing = bearing,
hab_p = hab_p,
mem_p = mem_p))
}Setup starting locations
# for the starting points of each individual buffalo
start_buffalo <- matrix(c(start_buffalo_df$x - xmin,
start_buffalo_df$y - ymin),
ncol = 2, nrow = nrow(start_buffalo_df))
# for random starting locations
start_unif <- cbind(runif(1e3,
min = ext(naive_predictions_cropped[[1]])[1],
max = ext(naive_predictions_cropped[[1]])[2]),
runif(1e3,
min = ext(naive_predictions_cropped[[1]])[3],
max = ext(naive_predictions_cropped[[1]])[4]))Simulating trajectories
n_indvs <- 3 # number of individual trajectories per starting location
n_steps <- 1000 # number of steps in the trajectory
n_ch <- 50 # number of proposed steps to be selected from
coef <- hourly_coefs # hourly coefficients
resc <- naive_predictions_cropped # resource selection (naive) raster stack
# most of the memory has declined to 0 by 500 locations
memory_period <- 500
memory_delay <- 24
spatial_sd <- memory_params_ALLoptim$mean_kde_sd
exp_gamma <- memory_params_ALLoptim$exp_gamma_all_optim
for(k in 1:length(start_buffalo_df$id)) {
stps_mem_list <- vector(mode = "list", length = n_indvs)
tic()
for(j in 1:n_indvs) {
stps_mem_list[[j]] <- simulate_ssf_memory(
n_steps = n_steps,
n_ch = n_ch,
coef = hourly_coefs,
xy0 = start_buffalo[k,],
starting_time = ymd_hm("2018-07-01 01:00", tz = "Australia/Queensland"),
step_time_increment = 60,
resc = resc,
memory_period = memory_period,
memory_delay = memory_delay,
# the KDE bandwidth (spatial sd) and temporal decay (exp_gamma) parameters
# here are the same for all individuals, but they could be indexed over
# a table of individual-specific parameters
spatial_sd = spatial_sd,
exp_gamma = exp_gamma)
}
toc()
animals_mem <- data_frame(
id = paste0("id", start_buffalo_df$id[k], 1:n_indvs),
track = map(stps_mem_list, ~ amt::track(
x = .$x,
y = .$y,
# we set time column to start at - which should allow enough time to
# exclude the warm-up period (e.g. 500 locations)
t = ymd_hm("2018-07-01 00:00") + hours(1:nrow(.)),
hour = .$hour,
angle = .$angle,
bearing = .$bearing,
hab_p = .$hab_p,
mem_p = .$mem_p
)))
sim_data_mem <- unnest(animals_mem, cols = track)
#check that the time lines up to the hour
head(sim_data_mem)
# return the locations to the original scale
#(i.e. add the xmin and ymin back to the x and y coordinates)
sim_data_mem_truexy <- sim_data_mem %>% mutate(x_ = x_ + xmin, y_ = y_ + ymin)
print(ggplot() +
geom_raster(data = ndvi_df,
aes(x = x, y = y, fill = NDVI_discrete),
alpha = 0.5) +
geom_point(data = sim_data_mem_truexy,
aes(x = x_, y = y_, colour = id),
size = 0.75, alpha = 0.1) +
geom_path(data = sim_data_mem_truexy,
aes(x = x_, y = y_, colour = id),
alpha = 0.75) +
scale_fill_brewer(palette = "Greens") +
scale_color_viridis_d("Sim ID") +
coord_equal() +
theme_classic() +
theme(legend.position = "none"))
filename <- paste0("outputs/simulated trajectories/CLR TS 2pDaily GvM MemNat/id_",
start_buffalo_df$id[k],"_GvM_mem", memory_period, "_", memory_delay,
"_", n_indvs, "ind_", n_ch, "ch_", n_steps, "stp_",
Sys.Date(), ".csv")
# print(filename)
write_csv(sim_data_mem_truexy, filename)
}## 107.78 sec elapsed## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## ℹ Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.## 109.94 sec elapsed
## 101.12 sec elapsed
## 319.37 sec elapsed
## 108.22 sec elapsed
## 99.77 sec elapsed
## 97.97 sec elapsed
## 306.05 sec elapsed
## 104.11 sec elapsed
## 101.6 sec elapsed
## 100.89 sec elapsed
## 306.89 sec elapsed
## 106.12 sec elapsed
## 101.79 sec elapsed
## 100.48 sec elapsed
## 308.5 sec elapsed
## 105.39 sec elapsed
## 102.8 sec elapsed
## 104.61 sec elapsed
## 313.34 sec elapsed
## 104.87 sec elapsed
## 97.94 sec elapsed
## 93.55 sec elapsed
## 296.89 sec elapsed
## 111.06 sec elapsed
## 102.08 sec elapsed
## 100.35 sec elapsed
## 313.81 sec elapsed
## 105.73 sec elapsed
## 98.21 sec elapsed
## 100.69 sec elapsed
## 304.83 sec elapsed
## 142.28 sec elapsed
## 98.95 sec elapsed
## 95.02 sec elapsed
## 336.46 sec elapsed
## 101.52 sec elapsed
## 96.58 sec elapsed
## 98.2 sec elapsed
## 296.5 sec elapsed
## 105.87 sec elapsed
## 98.67 sec elapsed
## 102.54 sec elapsed
## 307.35 sec elapsed
## 162.12 sec elapsed
## 95.7 sec elapsed
## 97.86 sec elapsed
## 355.91 sec elapsed
## 105.59 sec elapsed
## 97.2 sec elapsed
## 98 sec elapsed
## 301.11 sec elapsed
References
Signer, Johannes, John Fieberg, and Tal Avgar. 2017. “Estimating
Utilization Distributions from Fitted Step‐selection Functions.”
Ecosphere 8 (4): e01771. https://doi.org/10.1002/ecs2.1771.
Session info
sessionInfo()## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_New Zealand.utf8 LC_CTYPE=English_New Zealand.utf8 LC_MONETARY=English_New Zealand.utf8
## [4] LC_NUMERIC=C LC_TIME=English_New Zealand.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] Rfast_2.0.7 RcppZiggurat_0.1.6 formatR_1.14 scales_1.2.1 glmmTMB_1.1.8
## [6] clogitL1_1.5 Rcpp_1.0.10 ecospat_3.5 TwoStepCLogit_1.2.5 survival_3.5-5
## [11] viridis_0.6.2 viridisLite_0.4.2 matrixStats_1.0.0 patchwork_1.1.2 ggpubr_0.6.0
## [16] adehabitatHR_0.4.21 adehabitatLT_0.3.27 CircStats_0.2-6 boot_1.3-28.1 MASS_7.3-59
## [21] adehabitatMA_0.3.16 ade4_1.7-22 sp_1.6-0 ks_1.14.0 beepr_1.3
## [26] tictoc_1.2 terra_1.7-23 amt_0.2.1.0 lubridate_1.9.2 forcats_1.0.0
## [31] stringr_1.5.0 dplyr_1.1.2 purrr_1.0.1 readr_2.1.4 tidyr_1.3.0
## [36] tibble_3.2.1 ggplot2_3.4.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] backports_1.4.1 Hmisc_5.0-1 systemfonts_1.0.4 plyr_1.8.8 TMB_1.9.10
## [6] splines_4.2.1 TH.data_1.1-2 digest_0.6.31 foreach_1.5.2 htmltools_0.5.5
## [11] earth_5.3.2 fansi_1.0.4 magrittr_2.0.3 checkmate_2.1.0 cluster_2.1.4
## [16] tzdb_0.3.0 vroom_1.6.1 sandwich_3.0-2 timechange_0.2.0 colorspace_2.1-0
## [21] textshaping_0.3.6 rbibutils_2.2.13 xfun_0.39 crayon_1.5.2 jsonlite_1.8.4
## [26] lme4_1.1-32 zoo_1.8-12 iterators_1.0.14 ape_5.7-1 glue_1.6.2
## [31] PresenceAbsence_1.1.11 gtable_0.3.3 emmeans_1.8.5 car_3.1-2 abind_1.4-5
## [36] mvtnorm_1.1-3 DBI_1.1.3 rstatix_0.7.2 isoband_0.2.7 plotrix_3.8-2
## [41] xtable_1.8-4 htmlTable_2.4.1 units_0.8-1 foreign_0.8-84 bit_4.0.5
## [46] proxy_0.4-27 mclust_6.0.0 Formula_1.2-5 htmlwidgets_1.6.2 RColorBrewer_1.1-3
## [51] nabor_0.5.0 pkgconfig_2.0.3 reshape_0.8.9 farver_2.1.1 nnet_7.3-18
## [56] sass_0.4.5 utf8_1.2.3 tidyselect_1.2.0 labeling_0.4.2 rlang_1.1.0
## [61] reshape2_1.4.4 munsell_0.5.0 TeachingDemos_2.12 tools_4.2.1 cachem_1.0.7
## [66] cli_3.6.1 generics_0.1.3 audio_0.1-10 broom_1.0.4 evaluate_0.20
## [71] fastmap_1.1.1 yaml_2.3.7 ragg_1.2.5 maxnet_0.1.4 knitr_1.42
## [76] bit64_4.0.5 fitdistrplus_1.1-8 randomForest_4.7-1.1 nlme_3.1-162 pracma_2.4.2
## [81] biomod2_4.2-2 compiler_4.2.1 rstudioapi_0.14 e1071_1.7-13 ggsignif_0.6.4
## [86] bslib_0.4.2 stringi_1.7.12 plotmo_3.6.2 highr_0.10 lattice_0.21-8
## [91] poibin_1.5 circular_0.4-95 Matrix_1.6-5 nloptr_2.0.3 classInt_0.4-9
## [96] permute_0.9-7 vegan_2.6-4 gbm_2.1.8.1 vctrs_0.6.2 pillar_1.9.0
## [101] lifecycle_1.0.3 Rdpack_2.4 jquerylib_0.1.4 estimability_1.4.1 cowplot_1.1.1
## [106] data.table_1.14.8 raster_3.6-20 R6_2.5.1 KernSmooth_2.23-20 gridExtra_2.3
## [111] codetools_0.2-19 gtools_3.9.4 withr_2.5.0 multcomp_1.4-23 mgcv_1.8-42
## [116] parallel_4.2.1 hms_1.1.3 grid_4.2.1 rpart_4.1.19 coda_0.19-4
## [121] minqa_1.2.5 class_7.3-21 rmarkdown_2.21 carData_3.0-5 mda_0.5-3
## [126] sf_1.0-12 pROC_1.18.0 numDeriv_2016.8-1.1 base64enc_0.1-3 tinytex_0.48